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The response to a localized force provides a sensitive test for different models of stress transmission 
in granular solids. The elasto-plastic models traditionally used by engineers have been challenged by 
theoretical and experimental results which suggest a wave-like (hyperbolic) propagation of the stress, 
as opposed to the elliptic equations of static elasticity. Numerical simulations of two-dimensional 
granular systems subject to a localized external force are employed to examine the nature of stress 
transmission in these systems as a function of the magnitude of the applied force, the frictional 
parameters and the disorder (polydispersity). The results indicate that in large systems (typically 
considered by engineers), the response is close to that predicted by isotropic elasticity whereas the 
response of small systems (or when sufficiently large forces are applied) is strongly anisotropic. In 
the latter case the applied force induces changes in the contact network accompanied by frictional 
sliding. The larger the coefficient of static friction, the more extended is the range of forces for which 
the response is elastic and the smaller the anisotropy. Increasing the degree of polydispersity (for 
the range studied, up to 25%) decreases the range of elastic response. This article is an extension 
of a previously published letter fTl . 
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I. INTRODUCTION 

The description of the static phase of granular matter, 
much Hke its dynamical phases, is of significant current 
interest. Some studies focus on constitutive relations in 
order to describe the stress distribution in these materi- 
als whereas other investigations probe their microscopic 
characteristics, such as the nature of the force distribu- 
tions and correlations. 

Static and quasi-static properties of granular materials 
are commonly modeled by engineers using elasto-plastic 
[21 Ej and hypoplastic "4", "5] models. Other types of mod- 
els have been put forward (mostly) by physicists. Some 
of these models comprise hyperbolic partial differential 
equations for the stress transmission through a granular 
material O H [HI El [lOl E] . This is in contrast with the 
elliptic, non-propagating nature of the classical equations 
of static elasticity. It seems that this dichotomy in the 
modeling of granular statics started [3 [T2] in the context 
of the interpretation of experiments performed on conical 
piles [T5] [2] , where a pressure dip below the apex of the 
pile has been observed (the presence or absence of such 
a dip in conical piles and in wedges was found to depend 
on the construction method |15|). 

In order to obtain insights into this problem it is con- 
venient to consider the simpler geometry of a granular 
rectangular layer (or slab) resting on a solid support [16j , 
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not the least since this class of systems is well researched 
experimentally [T71 [TSl [H [201 [211 [21 [13] . Much hke in 
the case of piles, the cited and other experiments seem to 
be sending mixed messages concerning the 'correct' de- 
scription of stress transmission in granular slabs, as some 
render support to the elliptic descriptions whereas oth- 
ers are compatible with hyperbolic models. When the 
response to the application of an external force is linear 
in the force (for sufficiently small forces), the problem is 
tantamount to the study of the Green's function of the 
system [22j . 

The present article is devoted to a detailed study of 
the response of two-dimensional vertical slabs compris- 
ing polydisperse frictional disks subject to gravity and a 
localized, compressive external load applied at the top 
of the slab; it expands upon the results of [IJ in both 
content and detail. 

The dependence of the response on the magnitude of 
the external force as well as the interaction parameters 
and degree of disorder (or polydispersity) is investigated. 
The results suggest a resolution of the above mentioned 
controversy as well as the seemingly mutually incompat- 
ible experimental results. First, noticing that strongly 
anisotropic elasticity can exhibit hyperbolic-like features 
(see more below) it is convenient to consider the response 
within the framework of elasticity. Note that it is not 
claimed that the process leading to a given state is elas- 
tic (often it is not) but rather that the excess stress field 
induced by an external force in a given ( "prestressed" ) 
static state can be described by the equations of elas- 
ticity and often even by linear isotropic elasticity (for 
an isotropic reference configuration), even in the pres- 



ence of friction (see also the recent elastic model re- 
viewed in [21] )• More specifically, large systems subject 
to "small" applied forces (which is likely to be the case 
in many engineering applications), can be described by 
isotropic elasticity except in the near vicinity of the point 
of application of the force (where one expects strong in- 
duced anisotropy and irreversible rearrangements). Suf- 
ficiently far from the point of application the external 
force leads to practically no sliding or rearrangement of 
the grains. The latter events are rarer the larger the 
friction and therefore, as shown below, the linearity and 
isotropy of the response are enhanced by friction. 



Nonlinear/irreversible effects such as changes in the 
contact network and frictional sliding can naturally give 
rise to large anisotropy, yielding hyperbolic-like response 
(as already suggested in [HIS^). In relatively small sys- 
tems subject to sufficiently large forces, as in some ex- 
periments [HI [ISl HH HH [15], the size of the nonlinear 
and anisotropic domain induced by the external forcing 
can be comparable to the system's size, hence the stress 
transmission in such systems can be well described by hy- 
perbolic or near-hyperbolic (see below) equations. Simi- 
lar qualitative behavior is observed in both ordered (lat- 
tice) and disordered (polydisperse) configurations, the 
main difference being the range of magnitudes of external 
forces for which one observes elastic response. 



Note that qualitatively similar results were obtained 
in a study of the displacement response rather than 
stress [26j . using a similar numerical simulation. This 
study focused on the effect of the mean coordination 
number, controlled by the particle stiffness, rather than 
the effects of the frictional properties and the disorder 
which are the focus of the present work. Interestingly, an 
approach using a force network ensemble [.27j also seems 
to yield qualitatively similar results concerning the effect 
of friction, although a full physical interpretation of this 
approach is still lacking. 



The paper is organized as follows. Section II presents 
an outline of the elliptic (elastic) and hyperbolic de- 
scriptions of static granular matter and a comparison 
between the two. Results of some relevant experiments 
are reviewed in the same context. Sec. |III| provides a 
description of the simulation method. Sec. |IV| reports 
the results obtained for frictionless systems, and the 
subsequent sections describe more realistic simulations 
incorporating friction (Sec.|V| and disorder (polydisperse 
systems. Sec. 



VI I. The results of the simulations and 



their analysis enable the interpretation of the findings 
in different experiments: this comprises the content of 
Sec. |VII[ Sec. VIII offers a brief summary of the main 
results. 



II. ELLIPTIC VS. HYPERBOLIC 
DESCRIPTIONS 

A. Theoretical description 

As mentioned, there are two main classes of models 
that have been proposed for the description of the re- 
sponse of granular assemblies. Broadly known as the "el- 
liptic" and "hyperbolic" descriptions, their predictions 
differ both qualitatively and quantitatively. These re- 
spective predictions for the case of the response (excess 
pressure on the floor due to the external force, i.e., the 
pressure due to gravity in the otherwise unforced sys- 
tem is subtracted) to the application of an external ver- 
tical force at the center top of a slab are schematically 
illustrated in Fig. [l] The elliptic case is represented 
here by the results of isotropic elasticity [2E,- In this 
case, a single peaked response on the floor is expected 
[Fig. |l(a)| , its width being proportional to the depth 
(or height) of the slab. The shape of the peak is de- 
termined by the equations of elasticity, and depends on 
the boundary conditions at the floor [31 HOI US] (e.g., 
its rigidity and roughness). In contrast, hyperbolic de- 
scriptions [61 El [g [91 Uni dl] dictate that the stress 
propagates along characteristic directions, and two peaks 
are expected in 2D (a ring in 3D), their widths deter- 
mined by diffusive broadening due to disorder |30j and 
(consequently) proportional to the square root of the 
depth of the slab. This type of description has been 
shown [3 [HI [SI Uni E] to apply to (frictionless) isostatic 
systems [3T1 [321 [33] , i.e., systems in which Newton's first 
and third laws (and the boundary conditions) provide a 
number of equations that equals the number of unknowns 
(force components) thereby enabling the determination 
of the forces from these laws; it turns out these conditions 
determine an average coordination number for the parti- 
cle contacts and that the latter is the minimum required 
to maintain the stability of the packing. This notion 
cannot be directly extended to hyperstatic systems (i.e., 
with a higher coordination number). The isostatic case is 
a marginally stable configuration, which may correspond 
to a plastic material which is everywhere at incipient 
yield (for which classical plasticity results in hyperbolic 
equations [3]). We also note that recent models based 
on the description of the stress in terms of propagating 
force chains which may split and merge [311 133 1311 HZ] 
have been shown to correspond to elastic-like equations 
at large scales [33J [351 I3S]; a physical interpretation of 
the fields satisfying these equations is still lacking. 

It is important to note that, as already mentioned 
in [311251121], elastic systems can exhibit hyperbolic- like 
behavior. In |29j . we used a simple 2D model in order to 
demonstrate this property, i.e., a triangular lattice of har- 
monic springs in which the spring constant for the hor- 
izontal springs, Ki, is different from that of the oblique 
ones, K2. This model corresponds (in the continuum 
limit) to an anisotropic elastic system (the response of an 
anisotropic elastic infinite half-plane (2D) with uniaxial 




(a) "Elliptic" 



(b) "Hyperbolic" 



FIG. 1: Predictions of different models for the response of a 
granular slab to a localized vertical force applied to its top; 
pressure distribution on the floor supporting the slab. 



symmetry has been analyzed in detail in [36' ) . The qual- 
itative nature of the response depends on the degree of 
anisotropy: it is single peaked for values of K2/K1 near 1 
{K2/K1 = 1 corresponds to an isotropic system). The re- 
sponse is narrower than the isotropic one for K2/K1 < 1 
and wider for K2/K1 > 1. For sufficiently large val- 
ues of K2/K1, the response becomes double peaked, i.e., 
in qualitative agreement with the predictions of the hy- 
perbolic models. However, the equations of anisotropic 
elasticity are elliptic, except in the extreme anisotropic 
limit K2/K1 — > cx) in which the equations do become 
hyperbolic. This is consistent with the fact that the 
limit Ki —f corresponds to the absence of horizontal 
springs, and is therefore isostatic. This limit is simi- 
lar to the case of a square (or cubic) lattice of springs 
which (to linear order in the displacement) has no shear 
rigidity along the lattice directions [3^. Although this 
anisotropic model seems quite artificial, we show below 
that anisotropy arises quite naturally in more realistic 
descriptions of granular materials. 



B. A Review of Experimental Results 

This Section reviews some of the experiments in which 
the response of a granular slab to a localized force was 
measured. 

In the experiments described in [T8l|22], a localized ver- 
tical force was applied to a two-dimensional (2D) packing 
of photoelastic particles. Recall that stressed photoelas- 
tic solids viewed through cross-polarizers exhibit fringes 
whose density is proportional to the difference between 
the two principal values of the internal particle stress |39j . 
In the described experiment the applied forces had to be 
sufficiently large in order to observe the fringes; in prac- 
tice forces that were about 150 times the particle weight 
were used. The stress is concentrated at the interpar- 
ticle contacts. Post-processing techniques were used to 
extract the magnitudes of the interparticle forces from 
the experimental images [THl [551 ISH]- Force chains were 
prominent in all studied realizations. 

Three different types of packings were studied in [18^ 



155] . The most ordered was that of a triangular lattice 
of nominally monodisperse (equal) disks; less order was 
found in a packing of bidisperse disks (two different sizes), 
and a packing of pentagonally shaped particles was the 
most disordered of the three. The images from 50 con- 
figurations were averaged [T8j ]22\ for each type of pack- 
ing. The force profiles as a function of the horizontal 
(orthogonal to gravity) coordinate were plotted for dif- 
ferent depths in each of the slabs [TSJ [55]. These pro- 
files exhibited strong dependence on the particle shapes: 
the ordered lattice of monodisperse disks exhibited two 
prominent "force chains" along the lattice directions (60° 
with respect to the horizontal), a result which appears 
to be consistent with a hyperbolic model of propagat- 
ing forces (note, however, that a vertical force chain can 
be observed as well, and it is not anticipated by these 
models). As the disorder is increased, the force chains 
"fade out" , and for the random configuration of pentag- 
onal particles in 2D the measured force profile possessed 
a single peak, the width of which varied linearly in the 
depth; the latter result is compatible with the predictions 
of elasticity. It was proposed that granular materials may 
experience a crossover from a hyperbolic to an elliptic be- 
havior as the degree of disorder increases [18J. 

In another experiment, in which a rather small system 
was studied (of about 10 diameters in height), the slab 
comprised round-edged rectangular 2D particles |17j . In 
this case, the response appeared to exhibit a parabolic 
envelope. This type of force profile was predicted on the 
basis of models that assume an uncorrelated, diffusive 
propagation of the forces, such as the q-model [lOJ [S] 
(which was originally proposed to describe the statistics 
of force fluctuations). However, as mentioned, this sys- 
tem is quite small in terms of the number of constituents; 
furthermore, even in this case, it appears that the enve- 
lope may be narrowing down near the bottom of the slab. 
We shall return to the issue of system size below. 

In the simulations reported in [23j the displacement 
(rather than the force) response of a 2D packing of disks 
of three different diameters was measured (in response 
to a small displacement of a particle at the bottom of 
the packing). The displacement response is expected to 
be qualitatively similar to the force response (in isostatic 
systems the two are equivalent [31] [35]). The measured 
response (averaged over several hundred configurations) 
is single peaked, and the width of the peak appears to 
grow as the square root of the height for small system 
heights, and cross over to a linear dependence at larger 
heights (the total height of the packing was about 10 
diameters, which is quite small). The response function 
(scaled by the width) is fit quite well by a Gaussian. The 
same paper [23^ presents simulations of isostatic systems 
which show, as expected, a double peaked response. 

In the experiments described in ^^ [5D] , an external lo- 
calized force was applied at the top of 3D disordered slabs 
of sand (and other materials), and the response force pro- 
files (i.e., the actual forces minus the corresponding forces 
without the external load) on the floor were measured. 



In this case the apphed force was rather small (a few par- 
ticle weights), and the response was verified to be linear 
in the applied force. The measured response was single 
peaked, the profile width being proportional to the slab 
depth, as expected for an elastic system. However, the 
shape of the response was found to depend on the prepa- 
ration method: it was wider for loose packings (obtained 
by pulling a sieve through the packing) than for dense 
packings (obtained by filling the container layer by layer, 
pressing the packing after each layer is deposited). The 
response could not be fitted by isotropic elastic solutions 
for finite slabs PU". 

Unlike the highly disordered systems used in [T51 I22j . 
the experiments reported in [^T] concerned ordered 3D 
close-packed systems of spherical particles, arranged on 
FCC (Face-Centered Cubic) and HCP (Hexagonal Close- 
Packed) lattices |12]. The forces on the floor were mea- 
sured using pressure marks on carbon paper |21j : the 
applied forces used in these experiments were quite large 
(a few thousand times the particle weight). The results 
were averaged over several packings of nominally equal 
spheres. For shallow systems, three distinct peaks were 
observed in the case of the FCC lattice, and a sharp 
ring for the HCP lattice. This is consistent with a de- 
scription based on propagating forces |21j (which would 
correspond to a hyperbolic continuum description) . How- 
ever, for larger systems the peaks were considerably less 
distinct, the response being 'smoother'. This further in- 
dicates that a hyperbolic-like response may apply only 
to relatively small systems. For disordered (amorphous) 
packings [STJ , a single peak was obtained (a force impulse 
was used, since a large persistent force resulted in major 
rearrangements in the packing). 

In summary, qualitatively different types of response to 
an external force were observed in experiments on gran- 
ular assemblies, some of which seem to render support to 
elliptic models, others to hyperbolic or parabolic models 
of stress transmission. The degree of disorder seems to 
play an important role, as do the size of the system and 
the magnitude of the apphed force. The observed differ- 
ences also suggest that other parameters, such as the co- 
efficient of friction, may be important. In order to clarify 
the role of these parameters and the associated physical 
mechanisms, we performed extensive simulations of 2D 
granular slabs, as described below. 



III. SIMULATION METHOD 

In order to study the response of granular packings and 
its dependence on the particle properties in more detail, 
we employ the discrete element method [43l |44] (DEM) , 
a particle-based simulation method for granular materi- 
als also known as the molecular dynamics (MD) simu- 
lation method [45l |46l 113 |48]. This method has been 
employed in studies of atomic and molecular assemblies 
as well as granular systems, see, e.g., |35J [SHI [SI]) for 
recent reviews. In the present work, DEM simulations 



have been employed to study 2D systems composed of 
disks. The equations of motion for each particle are inte- 
grated using the Beeman algorithm |46j . which provides 
more accurate velocities than the commonly used Verlet 
algorithm [JB] . We verified that the use of a Gear 5- value 
predictor-corrector algorithm [48 does not result in any 
significant changes in the results. 



A. The Force Model 

Models for the interactions among solid grains are usu- 
ally based on contact mechanics [52l [53] . Following the 
work of Hertz (see, e.g., [2S]) and others, it is custom- 
ary to assume that the forces exerted by solid particles 
on each other are pairwise additive. Since for relatively 
rigid particles the typical deformation of a particle is a 
very small perturbation to its shape and size, it is com- 
mon to model these interactions as follows: the particles 
are fixed in shape and can overlap and the correspond- 
ing normal forces (perpendicular to the plane of contact 
among particles) are taken to be functions of the degree 
of this imaginary overlap. The tangential interparticle 
forces are taken to depend on the particles' rotations as 
well. 

Most natural or industrial granular materials comprise 
non-spherical particles. However, the latter are difficult 
to treat theoretically, and only few models for simu- 
lating their behavior have been suggested and studied 
(e.g., [54, 55, 56, 57l[58l[5g). Therefore the present pa- 
per specializes to disk shaped particles in two dimensions. 
The force model we use is essentially the same as that em- 
ployed in |3S] (for some other examples of force schemes 
commonly used in simulations, see, e.g., [60 j IHTl [62] ) . 

For spherical or disk-shaped particles, the overlap of 
two particles is defined by: 
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where Ri, Rj are the radii of the particles, r^ is the posi- 
tion of the center of mass of particle, i, and r^ = r^ — r-,-. 
In noncohesive granular materials, the particles are as- 
sumed to interact only when they overlap, i.e., when 
^ij > 0. For two frictionless elastic spheres, a classi- 
cal result by Hertz (see, e.g., [2S]) is that the force is 

3/2 

proportional to ^^ ■ , while for cylinders, it is linear in 
5y (up to a logarithmic correction, see, e.g., [63l [64]). 
In the simulations whose results are presented below the 
normal component (parallel to Vij) of the force is taken 
to be linear in the overlap (linear spring); this choice is 
not due to the fact that we consider cylindrical particles, 
but rather in order to simplify the application of theo- 
retical considerations. It may be further justified by the 
fact that we consider small deformations, in which case 
one may linearize the force-displacement law around a 
reference configuration |65j . 

Even for frictionless particles there can be dissipa- 
tion, e.g., due to viscoelastic normal forces that depend 



on the normal relative velocity fy (see, e.g., [MJ [55]). 
While hysteretic, rate independent dissipation models 
(see, e.g., [521 EZl EH]) may be more realistic than vis- 
coelastic models, their implementation is more compli- 
cated. Since we focus on the static response and small 
deformations the results should not be sensitive to the 
precise choice of dissipation mechanism. We therefore 
use a damping term which is linear in the relative ve- 
locity of interacting particles (linear dashpot). All in all 
the normal interparticle force exerted by particle j on 
particle i, f/^, can be expressed as: 
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where Vij is the unit vector in the direction of r^, v^. 
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Vi, • Vij is the relative normal velocity of 



the particles, v^ is a fixed damping constant, and H{x) 
is the Heaviside function, 
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Next, consider the friction-induced tangential forces. 
A simple model for friction [33] includes tangential 
springs for modeling static friction, often with velocity- 
dependent damping to facilitate relaxation to a static 
configuration. The springs are generated at zero length, 
which is also their rest length, when two particles come 
into contact and "snap" when the resulting frictional 
force, /-'", satisfies /^ > nsf^ , i.e., when the Coulomb 
limit is exceeded. Once a tangential spring is severed, the 
corresponding particles are allowed to slip with respect to 
each other and experience dynamic friction, /-^ = fijjf^. 
Note that the description of static friction in terms of 
a tangential spring is physically reasonable, and consis- 
tent with both contact mechanics on the macroscopic 
scale [52j [53] and microscopic approaches to the descrip- 
tion of friction (see, e.g., [SS] for a list of references on fric- 
tion) . The tangential springs arc of course quite stiff, giv- 
ing rise to very small, but measurable [70j displacement 
prior to slip. Tangential springs have often been used to 
model interparticle interactions in mean-field derivations 
of granular elasticity [7TJ [75] |73J |71]. More realistic de- 
scriptions of frictional contacts account for the history of 
contact deformations (see, e.g., [60 ] [61 ] [62 ] [67] ) . 

The above tangential spring-dashpot model in conjunc- 
tion with Coulomb's law of friction serve as the basic 
tangential force law in the present work, as in |43j . For 
simplicity, the coefficients of static and dynamic friction 
are taken equal each other: fj,D = fj,s = fJ. (the value of 
jiD does not directly affect the static configurations which 
are the main concern in the present work, although it may 
influence the evolution towards these states) . In order to 
employ a tangential spring-dashpot force model, a rel- 
ative tangential displacement has to be defined at each 
contact. To this end, note that the relative tangential 
velocity of two disks (in 2D) is given by: 



where Ui = Oi is the angular velocity and t^j = z x 
Yij. In the model used in this work the force exerted 
by a tangential spring depends on the relative tangential 
displacement, s^j, and is given by: 

where fc^ is the tangential spring constant and Sij obeys 
the following equation: 
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which ensures that the spring is always tangent to the 
plane of contact (a line in 2D), and does not exceed 
the length that corresponds to the Coulomb limit. In 
addition, s^j is set to zero if there is no overlap (i.e., 
^ij < 0). In practice, the tangential displacement is ob- 
tained by integrating the relative tangential velocity from 
the time the particles establish contact. When ^ifj > f/^, 
the length of the tangential spring is kept fixed in or- 
der to prevent the frictional force from exceeding the 
Coulomb limit. This is important in order to avoid un- 
realistic behavior when the contacts cease "sliding" (i.e., 
fT ^ ^fN^^ and revert to "sticking" (/^ < /x/^) [75]: if 
the integration of the tangential velocity is continued for 
a sliding contact, and the result used to further stretch 
the tangential spring, the result will be an unrealistically 
large force when the contact reverts to sticking. Note 
that the torques exerted on the particles determine their 
angular accelerations; the corresponding equations are 
part of the system of equations that are solved. 

The interactions of the particles with the walls are sim- 
ilar to the interparticle interactions. Since the walls con- 
sidered here are rigid and fixed straight lines (in 2D), the 
overlaps are calculated accordingly, with the same force 
models. The force constants for particle-wall interactions 
(kj/"", fc^^ii, lyjf"", v^^", ^^^") may be specified to be dif- 
ferent from those used for particle-particle interactions 
(which are taken to be the same for all particle pairs). 

In addition to interparticle forces and particle-wall 
forces, the force laws account for gravity. External forces 
and torques may also be applied to specific particles (see 
below) . 



B. Simulation Parameters 

The systems studied here consist of collections of poly- 
disperse disks, whose radii are uniformly distributed in 
the interval [R — S-R,R]. In order to facilitate the use of 
parameters from experiments performed on short cylin- 
ders, the cylinder thickness, W, and its volume density, 
p, are specified, and the particle masses are given by 
TTii — ttR^W p. We consider the case of homogeneous 
disks, for which the axial moment of inertia is given by 
li — 2"^i^f- The parameters used in this work cor- 
respond to experiments performed by the Duke group 
using 2D photoelastic disks [75": R — 3.75 • 10~^m; 
W = &.&■ 10-3 m; p = 1.15 • 10^ kg/m^. 



The simulation results are presented in terms of non- 
dimensional quantities defined as follows. The length 
unit is the mean particle radius R, the time unit is 
T = \l Kl g (where g is the magnitude of the accelera- 
tion of gravity) and the mass unit is the mean particle 
mass TO. In simulations which account for gravity, we use 
g = 9.8m/s2. 

The normal spring constant used for particle-particle 
interactions is fc^r = SOOOTTig/i?. This value is based 
on force-displacement measurements performed on the 
photoelastic disks mentioned above [§5]. Although the 
particles are cylindrical, measurements showed a force 
proportional to the | power of the displacement (pos- 
sibly indicating that the contact area is elliptic rather 
than rectangular) . For the typical deformations obtained 
in the experiments, a linear fit to the particle force law 
provides quite a good description |65j , which further jus- 
tifies our choice of modeling the normal force as that 
exerted by a linear spring with the above mentioned ef- 
fective spring constant. For particle- wall interactions, we 
use K^jf^ = 2fcjv. 

The damping coefficient is typically chosen to corre- 
spond to half the critical damping value for an individ- 
ual contact, Vc = 1\k (recall that fc and v are non- 
dimensional), i.e., we use vj^t — \/kM,T- This value 
was found to produce the fastest relaxation of the sys- 
tem towards static equilibrium and is irrelevant in the 
state of mechanical equilibrium itself. 

The time step used in the simulations was 

Af = 5- lO-^r = 5- lO-^yWs- We verified that 

decreasing the time step did not affect the results of the 
simulations. 

Other simulation parameters such as the ratio kx/kN 
(and kyp-^^ / k^J"^^ , which is taken to equal it), the coeffi- 
cients of friction /x, /i*^'' and the polydispersity 5 were 
varied in different runs, and their influence is discussed 
below. 



in time from zero to its desired value. The system is then 
relaxed again to static equilibrium. 

The main criterion used for testing whether the system 
has reached static equilibrium (at which time the simu- 
lation is stopped) is the kinetic energy per particle. We 
found that in order to obtain an accuracy of less than 1% 
for the forces on the floor in a given realization, which 
are used here to define the response of the system to a 
localized force, the system has to be relaxed to a kinetic 
energy of -E^. °^ = \Q~^^fhgR per particle, which is signif- 
icantly smaller than that used in previous studies (e.g., 
in [77 , the systems are relaxed to E'^ "'^ = 2 • lQ~^rhgR 
per particle). For polydisperse systems the response was 
coarse grained and averaged over a number of different 
configurations (Sec. VI I, since the forces exhibit strong 
fluctuations. In this case, sufficient accuracy was ob- 
tained with E^. "'^ = \G~^fhgR. We verified that several 
other criteria for static equilibrium (see |78j ) were satis- 
fied: the contact network was fixed for at least several 
hundred time steps (and each particle had at least two 
contacts), there were no sliding contacts, and the mean 
particle acceleration was less than lQ~^fhg. For friction- 
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(b)Forccs with a vertical applied force -Fext = brng at J. (gravity 
subtracted) . 



C. Initial Conditions and Simulation Procedure 

The initial configuration is produced by placing the 
particles on a triangular lattice with lattice constant 2R. 
The side walls and the fioor are placed at a distance R 
from the centers of the particles closest to the walls (so 
that for a monodisperse packing, Ri = R, these particles 
are tangent to the side walls and the fioor). The initial 
(translational and angular) velocities for all particles are 
set to zero. The system is then allowed to relax under 
gravity until it reaches static equilibrium, as described 
below. 

In order to study the response to a localized force, the 
resulting configuration (with all velocities, translational 
and angular, reset to zero) is used as an initial config- 
uration for a second run of the simulation, in which an 
additional external force (and/or an external torque in 
some simulations) is applied to a particle at the center of 
the top row of particles. The force is increased linearly 



FIG. 2: Forces obtained in a DEM simulation of a poly- 
disperse (5 — 0.25) frictional system with kr/kN ~ 0.8, 
fi — /i"''" — 0.2. The line widt hs a nd lengths are propor- 
tional to the force magnitudes, (a) After relaxation under 



gravity, (b) After relaxation with an additional vertical force 
-fcxt = 5fhg, where m is the mean particle mass, with the 



forces obtained in (a) subtracted. The lines drawn inside the 
particles indicate the rotation angles: in the initial configura- 
tion for the simulation (before the relaxation under gravity) 
the lines are vertical. An arrow denotes the position of the 
externally applied force. 

less systems {^ = /x^**" = 0), we found that reaching such 
small energies in a reasonable simulation time was very 
difficult due to the existence of slow global modes with 
a low dissipation rate, so that the following algorithm 
was used for accelerating the relaxation: the MD simula- 
tion was stopped at a higher energy (El °^ — lO^^fhgR). 
The configuration obtained at this stage was used as a 
reference state around which (since the interactions are 
harmonic in the static case) the linearized equations of 



motion were used iteratively to obtain a system which is 
closer to static equilibrium. In each iteration the equa- 
tions were solved by matrix inversion, and the connectiv- 
ity of the particles was updated to ensure that there was 
no tension. The iterations were stopped when the maxi- 
mal relative particle displacement was less than 10~^'^R 
(which is essentially the numerical accuracy). The typi- 
cal kinetic energy obtained when using the resulting con- 
figuration as an initial configuration for the full DEM 
simulation was less than Ek = lO^^^fhgR per particle. 
This configuration was then used for calculating the in- 
terparticle forces. 

An example of the interparticle forces obtained in a 
typical simulation run is presented in Fig. [2] In order 
to calculate the response to an applied force, the forces 
due to gravity alone (without the applied force) are sub- 
tracted (vectorially) at each contact. 



IV. RESPONSE OF FRICTIONLESS SLABS 

Consider first frictionless particles, interacting by uni- 
lateral ("one-sided") springs, which apply force only 
when compressed, modeling repulsive-only particle in- 
teractions 123], see Eq. (pi. As shown in [25], the ap- 
plication of a localized force at the top of the packing 
can lead to rearrangements in the contact network: hor- 
izontal springs in the region below the point of appli- 
cation of the force are opened (as also observed in [73] 
for a pile geometry). The force chains in this system 
are qualitatively similar to those obtained with regu- 
lar ("two-sided") springs. However, the force magni- 
tudes vs. the horizontal coordinate at different depths 
are expectedly in better agreement with experiment in 
the case of one-sided springs [181 122] than that of regu- 
lar springs |29j . since the former model is more realistic 
for granular systems. The stress distribution obtained 
for unilateral springs |29j is anisotropic, and shows two 
peaks at the floor. This anisotropy is obviously related 
to the existence of a region of open horizontal contacts 
where the anisotropy is large and the hyperbolic limit 
applies [29|. This anisotropy is not present in the sys- 
tem without the application of the external forces; it is 
induced (through the changes in the contact network) 
by the applied force [23]. The induced changes in the 
contact network obviously qualify as a nonlinear process. 

These changes in the contact network may be mod- 
eled as a nonlinear extension of the linear elastic con- 
tinuum behavior obtained in a network of harmonic 
springs. While the (possibly position-dependent) elas- 
tic moduli in linear elasticity are time-independent ma- 
terial properties, a possible extension is to introduce a 
stress history dependence of the elastic moduh (i.e., the 
anisotropy induced by the opening of contacts in certain 
regions may be considered a result of an attempted tensile 
stress in those regions). A similar type of stress-induced 
anisotropy has been suggested in the context of plastic 
models of soil mechanics 80 1 as well as in nonlinear elas- 



tic models [53]. If the particle positions do not change 
significantly, so that only the contact network is modified 
in response to the appHed stress, the behavior can possi- 
bly be modeled as "incrementally elastic" . Under certain 
conditions (corresponding to plastic yield) , the system is 
no longer able to support the applied stress without a 
major rearrangement of the particles. Incipient plastic 
yield may be related to a local extreme anisotropy typi- 
cal of a marginally stable isostatic configuration. 



A. Dependence on the Applied Force: Crossover 
from Hyperbolic-Like to Elliptic Response 

In order to examine the changes in the contact net- 
work in more detail, we performed DEM simulations, as 
described in Sec. Ill of a system similar to those dis- 
with different applied forces. We focused 

the effect 



cussed in [29] 

on systems of 15 layers of 60 particles each; 

of system size is discussed below. The force response on 

the floor (subtracting the effect of gravity) is shown in 

Fig. [3J A crossover from a single peaked to a double 

peaked response occurs as the applied force is increased. 

The changes in the contact network corresponding to 
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FIG. 3: The response of frictionless ordered systems for dif- 
ferent applied forces, Fcxt- 

the systems of Fig. |3| are shown in Fig. [4] For a suffi- 
ciently small force (not shown), the contact network is 
unchanged, and the response is fully elastic (see Sec. [V] 
for a discussion of the linearity of the response). As the 
force is increased, horizontal contacts are opened in a 
teardrop shaped region below the point of application 
of the force, whose size increases with the force and de- 
creases with friction. As mentioned above, in this re- 
gion the extreme anisotropic limit of elasticity applies. 
Whe n the "te ardrop" is sufficiently far from the floor 
[Fig. 4(a)||(b)| , the response at the floor is single peaked. 



as the changes induced by the force can be considered 
to be localized. Otherwise, the anisotropy induced by 
the external force reaches the floor (Fig. 4(c)| inducing 
a double peaked response (the crossover actually occurs 
at slightly smaller forces, at which the "teardrop" almost 



reaches the floor). 



(a)i^ext = l.bmg. 



IbJT^ 



: 3mg. 



(c)Fext = 6mg7 



FIG. 4: Changes in the contact network in an ordered fric- 
tionless system for different applied forces. The central third 
of the system is shown. Thick lines connecting particle cen- 
ters indicate contacts opened due to the application of the 
force; thin black lines represent contacts which are unaffected 
by the external force. 



B. Dependence on System Size 

The dependence of the crossover on the size of the sys- 
tem is important for the interpretation of experimental 
results. The depth and maximal width of the "teardrop" 
as a function of the applied force, for several system sizes, 
are plotted in Fig. |5] As the "process" of opening of 
contacts with increasing applied force is nonlinear (since 
the effective elastic properties of the system are modi- 
fied in the process), it is difficult to model it analytically. 
However, Fig. [5] clearly demonstrates that the size of the 
"teardrop" only slightly decreases as the depth of the slab 
increases. 

This indicates that for large systems (as typically found 
in nature and in engineering application), this change 
in the contact network gives a finite size correction to 
the response [T]. For identical external loads, very deep 
systems should exhibit a single peaked response whereas 
relatively shallow ones, as studied in some experiments 
(Sec. 



II B I, should exhibit two peaks. This is the main 



reason that different experiments, using differently sized 
samples, yield qualitatively different results. This obser- 
vation also applies when frictional interactions are ac- 
counted for, as described in Sec. IVl 



appreciable and the localized force would lead to signifi- 
cant rearrangements of the particles) . 

The independence on the particle stiffness may be un- 
derstood as follows: in the reference configuration, con- 
tacts are compressed due to gravity as well as the rigid 
walls and floor, i.e., the stress, and in particular its com- 
ponent parallel to the horizontal contacts, axx, is com- 
pressive (the same condition may of course be obtained 
by applying an external pressure to the system) . The ex- 
ternal force acts in the opposite direction, i.e., it attempts 
to give rise to tensile forces. Therefore, the opening of 
contacts is due to a "competition" between the compres- 
sion due to gravity and the tensile (excess) forces due to 
the applied force, and is therefore determined directly by 
the stress, rather than the strain (which does depend on 
the particle stiffness, of course), provided that the geom- 
etry is not signiflcantly affected by the applied force (as it 
would for particles which are sufficiently soft to allow sig- 
nificant overlap) . This justifies our choice of the particle 
weight as the unit of force (rather than a scale based on 
the particle stiffness and size). The same behavior may 
persist even in the limit fcjv ~* oo under the same bound- 
ary conditions provided that the particles in the reference 
state are still in contact along the horizontal direction. 
However, as the stiffness increases, this requirement may 
be hard to comply with, as the particles must fit exactly 
between the walls in order to remain in contact: in the 
limit fcjv —^ oo most horizontal contacts may be absent al- 
ready in the reference state, rendering it isostatic (there- 
fore, extremely anisotropic), and this would lead to a 
hyperbolic-like response even for an infinitesimal applied 
force, as suggested in, e.g., |71 |S1 IB UHl EI] • Therefore the 
above 'limit of infinite stiffness' should be understood as 
'large but finite' stiffness else isostaticity comes into play. 
Furthermore, one must bear in mind that frictionless par- 
ticles are a rather artificial idealization in the context of 
granular materials. Indeed, as we show in the next sec- 
tion, the presence of friction affects the behavior of the 
system very significantly and, in particular, renders it 
much less sensitive to the contact network. 



C. Dependence on Particle Stiffness 

It is clearly important to consider the effect of the stiff- 
ness of the particles. We performed simulations with 
different values of kjsi, and found that the size of the 
"teardrop" did increase with the stiffness of the parti- 
cles, but appeared to saturate for k^ ^ 2000mg/R. The 
results presented here pertain to k^ — SOOOmg / R. The 
response is also quite insensitive to the choice of stiff- 
ness for fcjv ^ lOOOmg / R: at a; = 0, it only changes 
by about 2% in the range lOOOmg/R < fcjv < lO^mg/R. 
We conclude that the crossover from a single peaked to 
a double peaked response is essentially independent of 
the particle rigidity (provided that the particles are suf- 
ficiently stiff; otherwise the particle "overlaps" may be 



V. EFFECTS OF FRICTION 

Dependence on the Applied Force: Friction 
Increases the Linear Range 



As mentioned in Sec. |IV| for sufficiently small external 
forces no changes are induced in the contact network, so 
that the response is expected to be linear in the applied 
force |81j . This is shown in more detail in Fig. m\ which 
presents the response at the floor (vertical force on the 
floor) at x = (below the point of application of the 
force) as a function of the applied force, in both friction- 
less and frictional systems (with kr/k^ — 0.8; the effect 
of this parameter is discussed below) . The results shown 
here were obtained with frictionless walls (/x'"'*" = 0; see 
Sec. Ill), since for weak applied forces, frictional walls can 
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FIG. 5: The depth and maximal width of the "teardrop" in units of the particle radius, R, in ordered (monodisperse) frictionless 
systems (Fig. p|, as a function of the applied force, normalized by the weight of a particle, for different system sizes. 
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FIG. 6: The response at a; = of ordered systems, normalized 
by the applied force, F^^t, vs. the applied force (given in units 
of the particle weight, mg), for different coefficients of friction 
fi, frictionless walls (/i™^" — 0) and kr/kN ~ 0.8. 



support some of the force, inducing nonlinearity in the 
response at the floor. A Unear (elastic) range is observed 
for a sufficiently small applied forces even for frictionless 
particles, due to the fact that the particles are slightly 
deformed by gravity, and a small force does not cause 
the contacts to open. 

As shown in Fig. |6] friction has a significant effect on 
the linearity of the response: The extent of the linear 
range is significantly larger (by almost an order of mag- 
nitude; notice the logarithmic scale on the horizontal axis 
in Fig. m\ in frictional systems, so that elasticity is en- 
hanced by friction. 

The effect of friction on the response is shown in more 
detail in Fig. [7] which presents the response profile on the 
fioor for different applied forces, for /i = 0.2 and /i = 1 
(compare to the frictionless case presented in Fig. Isl. 
Notice that the force for which the crossover from a sin- 



gle peaked to a double peaked response occurs rapidly 
increases with friction (this is also observed in Fig. l6|, 
so that friction renders the response closer to that ex- 
pected from isotropic elasticity (this is further discussed 
below). For /i = 1 no such crossover is observed even for 
the largest force shown, F^xt = 50mg (a different type of 
crossover occurs for fi — 1 for large forces, as described 
in Sec. VE). Note that much larger forces may induce 
major rearrangements (i.e., plastic ffow), while our fo- 
cus here is on the solid-like regime, in which the particle 
displacements are small. 



B. The Effect of Friction on the Contact Network 

To gain a better understanding of the effect of friction 
on the response, we examine the changes in the contact 
network. We first recall that such changes occur at sig- 
nificantly larger applied forces in the frictional case (by 
almost an order of magnitude): the first contact is opened 
when Fcxt = 0.75mg in the frictionless case, F^xt — '^rng 
for /i — 0.2, and i^oxt — Qrng for /i = 1. This is the 
origin of the extended linear range observed in frictional 
systems. 

The effect of the applied force on the changes in the 
contact network in a system with /i = 0.2 is shown in 
Fig. Is] (compare to Fig.|4|. For the same force, the region 
of open contacts is considerably smaller in the frictional 
case than in the frictionless case [compare Fig. 



Fig. 8(a) . In addition, this region reaches the ffoor only 



4(c) and 



for a much larger force [compare Fig. 4(c)|and Fig. 8(c) 



This may explain the fact that the crossover to a dou- 
ble peaked response occurs at larger forces in the fric- 
tional case, since this type of response is related to the 
anisotropy in this domain. 

The above observations of the changes in the contact 
network seem to explain the effect of friction on the range 
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FIG. 7: The response of ordered systems, with different coefficients of friction fi, to applied forces (-Foxt) of different magnitudes, 
frictionless walls (/i"''" = 0) and kr/kN = 0.8 (compare to Fig. Isk. 



l£^)i'cxt = 6rMy. 



(b)i'ext = lonig. (c)i'ext = SOrng. 



FIG. 8: Changes in the contact network in a frictional system 
with n — 0.2, with different applied forces Foxt, frictionless 
walls (^""'" = 0) and kr/kN = 0.8. The various lines are 
explained in the caption of Fig. |4] 



of linearity and magnitude of the crossover force. How- 
ever, the source of this effect, namely, the reason for 
the smaller changes in the contact network in the pres- 
ence of friction, still needs to be elucidated (see below). 
In addition, for a larger coefficient of friction, the ob- 
tained response does not seem to be compatible with the 
changes in the contact network: as mentioned above [see 
Fig. |7(b)| , for /i = 1 the response remains single peaked 
(and actually its shape becomes sharper) even for rather 
large forces (see however Sec. VE), while for the same 
force, the response for /x = 0.2 is double peaked. The 
changes in the contact network for /i = 1 are shown in 
Fig. [9) While the region of open contacts does reach the 
floor (even for a smaller force than in the case of /i = 0.2), 
the response remains single peaked. 
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FIG. 9: Changes in the contact network in a frictional system 
with /i = 1, with different applied forces Fext, frictionless walls 
^^waii ^ Q^ ^^^ kr/kN = 0.8. 

The results obtained for fi — 1 suggest that in the 
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FIG. 10: The response of ordered systems for a large coeffi- 
cient of friction {fj, = 10) and different values of the ratio of 
tangential to normal stiffness, kr/kN, with _Foxt = 15mg. 



presence of friction, the anisotropy of the region of open 
contacts is greatly reduced (compared to the frictionless 
case). When /i is sufficiently large (and Fcxt sufficiently 
small), its value is essentially immaterial, since sliding 
is never reached. In this case, the ratio of tangential to 
normal stiffness, kx/kj^ (see Sec. Ill A I, determines the 
response: Fig. [TO] presents the response obtained with 
-Fcxt = I5mg for /x = 10 (which is practically equiva- 
lent to /i = 00, as used in [26], since no sliding occurs), 
with different values of kx/kj^. A crossover from a sin- 
gle peaked to a double peaked response occurs with de- 
creasing kr/k^ (the crossover occurs at kr/k^f ~ 0.3). 
For larger values of kr/kN, the response is nearly inde- 
pendent of kx/k]^. This phenomenon is further studied 
immediately below. 
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FIG. 11: A spring model with normal and tangential springs. 



C. A Model for Large Friction: the Role of 
Tangential Stiffness 

In order to understand the dependence of the response 
on kx/kN, we consider a 2D spring network model sim- 
ilar to the one described in [291 (i.e., a triangular lattice 
with different spring constants for horizontal springs, /c", 
and oblique springs, fcj). In addition, this model (see 



Fig. 11 1 contains tangential springs (as used in our DEM 
simulations; Sec. Ill A), with spring constants which are 
different for horizontal, k[, and oblique contacts, /cj [SI]- 
Note that if one considers Hertzian interactions, the ex- 
istence of anisotropic prestress (e.g., due to gravity) may 
indeed result in different spring constants in the horizon- 
tal and oblique directions. This model does not incorpo- 
rate sliding, so that it corresponds to /x = oo. Similar 
models were studied in [TH |S2 US], and more recently 
in [36ll74[l84]. where equal spring constants were used in 
all lattice directions. 

The elastic moduli corresponding to the long- 
wavelength limit of this model can be calculated as fol- 
lows: to leading order in the relative particle displace- 
ments Uij (which is the relevant order for linear elastic- 
ity), the elastic energy of the system is given by: 



£;°^ = 



2 ^^ 



k7j [h 



+ K, [u,j 



u, 



,)f^: 



(5) 

where the sum is over nearest neighbors, and the values 
of A:" and fc* are taken according to Fig. fill and the 
above description. In order to obtain the elastic moduli 
corresponding to the continuum limit of the considered 
system, an affine deformation (which is appropriate for 
a lattice configuration) is defined by a symmetric, uni- 
form strain field: Uijaijc^t) = eafsrijfs- Note that when 
tangential forces are present, the stress is not necessarily 
symmetric (as assumed in classical elasticity); however, 
the micropolar terms due to this asymmetry are of higher 
order in the strain. We verified (in the DEM simulations) 
that they are very small, even near the point of applica- 
tion of the external force, where the magnitude of the 



antisymmetric part of the stress is only a few percent of 
the pressure. 

Using the notation of [55] : 



E'' 




(6) 



where the superscript, T, denotes the transpose, one ob- 
tains the following elastic moduli: 



« = ^ (8fcr + ^2 + 3fc*) 

d2 
C — ITT (3^2 ~ 3fc2J 

d = ^{Qk^+^kl+2kl) 



2A 



(7) 

(8) 

(9) 

(10) 



where A — 2^/3R^ is the area of the unit cell. These re- 
sults are consistent with those obtained in [36 for normal 
springs only, but different from the model with bending 
interactions introduced in |36j . 

In a domain of open horizontal contac ts o ne has fc" = 
k\ — 0. It was already shown in Sec. IV that in the 
absence of tangential forces (fc| = 0), this system corre- 
sponds to the extreme anisotropic limit. However, since 
the oblique tangential springs apply forces which have 
horizontal components, they can (at least partially) com- 
pensate for the absence of normal horizontal springs, and 
therefore significantly decrease the anisotropy. 

Otto et al. [33] present continuum elastic solutions for 
an anisotropic infinite half plane subject to a localized 
force. They found a criterion for the elastic moduli at 
which a crossover occurs from a single peaked response 
to a double peaked one ^: two peaks are expected 
for r = ^ [ab — c{d + c)] < 0. Note that this crite- 
rion refers to an infinite half-plane rather than a slab 
of finite width, which is more appropriate for describ- 
ing our simulations. In addition, this model is homoge- 
neous, while the region of open contacts (even when it 
reaches the floor) is roughly a triangle below the point 
of application of the force. Nevertheless, the estimate 
obtained using this criterion fits our result for the finite 
slab quite well: for the model used here (with no hori- 



zontal springs). 



_ 13^ + 1013-3 
f3^+6P+9 



where /3 = k\/k2- Hence 



(in the physically relevant range /3 > 0), two peaks are 
expected for (3 < 0.2915, which is consistent with the 
value kx/k^ — 0.3 obtained in the simulation (Fig. 10 1. 



For an infinitesimal tangential load applied to a sys- 
tem composed of elastic spheres, the Cattaneo-Mindlin 
model [53 yields fcx/fcjv — 2-v ' "^^^^re v is the Pois- 
son ratio of the spheres. For the range of positive Poisson 
ratios {Q <v < 0.5; notice that this is a 3D Poisson's ra- 
tio), this implies that g < kx/k]^ < 1, so that the min- 
imum value of kx/kpf is well above the crossover from 
two peaks to one. In most of the simulations presented 
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FIG. 12: The response of ordered systems for kr/kN = 0.8 
and different values of the coefficient of friction, with Fcxt = 
15mg. 
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FIG. 13: The interparticle forces in ordered systems with dif- 
ferent coefficients of friction (p""'" = fi,kT/kN = 0.8), with 
-Foxt ~ 15mg. The effect of gravity has been subtracted. The 
central third of the system is shown. Line widths and lengths 
are proportional to the force magnitude. 
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in this work, we use kx/kj\j — 0.8, which corresponds to 
a realistic Poisson ratio of about 0.3. 



FIG. 14: The vertical stress a^z in the systems shown in 
Fig. |13| calculated using a Gaussian coarse graining function 
with a coarse graining width w = 2R. The effect of gravity 
has been subtracted. 



D. More on the effects of friction 

The above results show that static friction acts to re- 
tain an "effective connectivity" between the grains when 
the horizontal contacts are disconnected, and renders the 
system more isotropic than one would have naively antic- 
ipated. However, the frictional forces are limite d by the 
Coulomb condition (/-^ < fif^; see Sec. Ill A I, so that 
if /i is too small (or F^xt large), sliding occurs. This in- 
creases the anisotropy, since not all the tangential springs 
can exert forces as large as predicted by the above model 
(in which ^ = oo) and leads to a crossover to a double 
peaked response. This crossover is shown as a function 
of /i (for Fext = I5mg) in Fig. 12 



In considering the effect of the coefficient of static fric- 
tion, note that unlike the coefficient of dynamic friction, 
which is typically smaller than 1, the effective coefficient 
of static friction (which determines the onset of sliding) 
may be significantly larger than 1 even for rough spherical 
particles, and certainly for irregularly shaped particles 
(which can interlock and prevent their relative rotation) . 

The effects of friction are also evident in the forces 
(Fig. [T3| as well as the vertical stress component azz 
(Fig. |14[ ), calculated using the expression presented 
in [29, i85j . In particular, the reduced anisotropy is ap- 
parent in the stress field (the stress field for /i = 1 is 
quite similar to that obtained in the case of an isotropic 
harmonic lattice [23] )• 

The most prominent effect of friction on the force 
chains (Fig. 13 1 is that the forces (and hence the chains) 
no longer need to be aligned with the (here, triangular) 
lattice direction, as they are in the absence of tangential 
forces. In particular, a vertical force chain is evident for 
/x = 1. Such a chain is actually observed in the exper- 



iment [m m] . The inclusion of friction does provide a 
better fit to the experimental results obtained in [18l [22] 
for the monodisperse ordered packing: the force magni- 
tudes vs. the horizontal coordinate at different depths in 
DEM simulations with frictionless and frictional (/z = 1; 
for the particles used in the experiments, ^ ~ 0.94 [76] ) 
are shown in Fig. [15] We have also been able to re- 
produce [HS] force chains along non-lattice directions ob- 
served in experiments with an applied force at oblique 
angles 22J, using similar simulations of frictional parti- 
cles (with some polydispersity; see Sec. |VI| below). 



E. Symmetry Breaking for /i = 1 

As described above, for a relatively small coefficient 
of friction we observe a gradual crossover from a sin- 
gle peaked to a double peaked response as a function 
of the magnitude of the externally applied force, much 
like in frictionless systems, the main difference being that 
it occurs at larger values of the applied force [compare 
Figs. [3] and [7(a)] , although the shape of the response is 
not identical. We reiterate that for the frictional case, 
this crossover can be explained as follows: in the limit 
^ — > 00, the induced anisotropy due the changes in the 
contact network (the opening of contacts; see Figs.[4[and 
^ is compensated by the tangential forces, which renders 
the response single peaked for kx/kj^ > 0.3. For finite fi, 
this compensation is limited by the Coulomb condition, a 
limitation which becomes more significant as the applied 
force is increased: for /i = 0.2, the region of open con- 
tacts (the "teardrop" ) reaches the floor for F^xt — 20mg, 
which in this case seems to be the crossover force [see 
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FIG. 15: The norms of the interparticle forces, |f|, vs. the horizontal position, x, in ordered systems, with applied force 
Fcxt ~ 150mg. The legend indicates the depth measu red from point of application of the force, in particle radii (compare 



to the experimental measurements shown in 
kr/kN = 0.8. 



(a) Frictionless particles, (b) Frictional particles with /i 
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Fig. 16(a)| . Sliding does occur in the system (mainly near 
the point of application) even for smaller forces, and slid- 
ing contacts spread gradually downward with increasing 
Fcxt (possibly due to the fact that the displacements in- 
duced by the applied force decay with the distance from 
it). It therefore appears that the crossover for fi = 0.2 
is associated with the gradual reduction of the compen- 
sating effect of the tangential forces due to sliding. For 
larger forces (i^oxt ^ 40TOg), a small third peak appear 
at X = (below the point of application) , which appears 
to be related to a reorganization of the sliding contacts 
within the "teardrop" . 

For fi = I, however, we observe quite a different be- 
havior. As described in Sec. |V A[ the response remains 
single peaked (and even becomes sharper) for rather large 
forces, much beyond those for which the crossover is ob- 
tained for /J, = 0.2 [see Fig. |7(b)| . However, when the 
force exceeds i^ext = 91mg, we obtain a sharp transition 
to a very different response shape [see Fig. 16(b)| , which 
is very asymmetric (unlike the symmetric double peaked 
response obtained for lower friction at large forces). The 
sign of symmetry breaking appears to be related to the 
transient vibrations which are excited by the application 
of the force; for some very close values of the force [e.g., 
Fext = Q^mg, shown in Fig. |16(b)| we obtain almost ex- 
actly the same rescaled response reflected with respect to 
the vertical axis through x ~ 0. The shape remains qual- 
itatively similar for much larger Fcxt- The transition to 
an asymmetric response is accompanied by a significant 
change in the contact network (which becomes asymmet- 
ric), and by a marked reduction in the number of sliding 
contacts. 

This instability at large /i appears to be dominated by 
sliding, rather than by changes in the contact network. 
Its source (and in particular the breaking of symmetry) is 



yet unclear, although it may be related to the well-known 
phenomenon of shear banding. It also bears some simi- 
larities to elastic buckling (however, it is clearly outside 
the linear elastic regime). 



F. A Pile Geometry 

While the work presented here is focused on the re- 
sponse of a granular slab, we also examined the effect 
of friction in a pile composed of 11 layers of monodis- 
perse disks arranged on a triangular lattice, prepared as 
described in Sec. IIIC (similar to the frictionless piles 
studied in [79]). These piles are quite unrealistic, since 
the sides of the pile are at 30° to the horizontal, which is 
larger than the realistic angle of repose for disks. There- 
fore, the edges of the pile have to be supported by side 
walls (the construction of stable piles on an uncorru- 
gated floor requires the introduction of rolling friction, 
or rolling resistance [57]). 

The forces and the corresponding vertical stress com- 
ponent a^z in the pile are shown in Figs. 17(a)||(d) As in 
the case of a slab (Figs. 13 and 14), the effect of friction 



is rather large. In particular, as shown in Figs 
the pile of frictionless disks exhibits a dip in the stress un 
der the apex, but the pile of frictional disks does not. The 
dip in the frictionless case is due to the open horizontal 
contacts in the central region of the pile; see Fig. 17(e) 



In the frictional case, the size of the region of open con- 
tacts in the central region of the pile is quite similar, 
but its shape is different [Figs. 17(e)|(f)] . The absence of 
the dip results from the reduced anisotropy due to the 
frictional forces, as discussed above for a slab geometry. 
We therefore conclude that the presence or absence of 
a dip in a pile geometry is determined by the degree of 



14 




ext 



0.15 

0.1 

0.05 



-0.05 

-0.1 

-0.15 

-60 -40 



^ 




- - Fex,=92 mg 
F^^ =95 mg 

- - Fex,=200 mg 



-20 
X/R 

{b)M = 1- 



20 



40 60 



FIG. 16: The response of ordered systems, with different coefficients of friction /i, to an apphed force (Fcxt) of different 
magnitudes, frictionless walls (/i™^" — 0) and kr/kN ~ 0.8. 



1 1 1 1 1 1 1 1 1 1 



(a)Forces, fj, = 0. 



.■./.'.'.".'VVVVVVVVVVVVVVVV. 

(b)Forces, fj. = 0.94. 



F' 



il 



w^ ^i 



c)Vertical stress, /i = 0. (d)Vertical stress, /i = 





(c) Contact network, /.i = 0. (f)Contact network, /i = 0.94. 



FIG. 17: (a)l(b) the interparticle forces in a pile of monodis- 
perse disks under gravity, with and without friction (in the 
frictional case, /i — 0.94, /x""" = 0.35 and kr/kN ~ 0.5). Line 
widths and lengths are proportional to the force magnitude; 
(c)l(d) the vertical stress ct^z in the same systems, calculated 
using a Gaussian coarse graining function with a coarse grain- 
ing width w = 2R; (e)|(f) the contact network: particles in 
contact are connected by lines. 



anisotropy in the mechanical properties of the pile (which 
is not always simply related to the anisotropy of the con- 
tact network, as shown here in the frictional case). These 
properties are expected to depend on the way the pile is 
prepared, as observed in [15]. 



G. Effects of an applied torque 

It is interesting to note that the force chains are also 
sensitive to additional applied torques. Fig. [18] shows 
the force chains obtained in systems of 25 x 13 slightly 
polydisperse particles {6 = 10"'^) with kx/kN = 0.5, 
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FIG. 18: Force chains in 2D packings of slightly polydisperse 
frictional particles {S = 10~^), with a force Fgxt = lOffig ap- 
plied at 30° to the horizontal, and different applied torques 
in the clockwise direction (indicated below each figure). The 
effect of gravity has been subtracted. The same realization 
of the packing was used in all cases. The region shown is the 
central third of the upper half of the system. 



30° to the horizontal, with different applied torques (as 
may occur if the line of action of an external force ap- 
plied to a particle does not pass through its center of 
mass). As shown, additional torques can have quite a 
large effect on the obtained force chains. Although the 
effect on the stress field (in particular, its asymmetric 
part) is typically confined to a region close to the point 
of application, this may indicate that for some effects it 
is important to extend the continuum description to in- 
corporate torques and rotations, i.e., use a micropolar or 
Cosserat continuum model [88 ] l89 ] l90]). 



VI. EFFECTS OF DISORDER 

While the study of ordered systems of monodisperse 
grains allows for simpler theoretical modeling, real gran- 
ular systems are never composed of perfectly identical 
particles. Even disks of nominally equal sizes (as used, 
e.g., in [T5J 122]) exhibit a certain polydispersity. Most 
real granular matter exhibits considerable polydispersity 
(as well as a variability in particle shapes). It is therefore 
important to verify that the results obtained for monodis- 
perse system, as described above, are not limited to this 
particular, idealized case. We therefore performed simu- 
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FIG. 19: The response of frictionless and frictional 
{fi = /i""'" = 0.2, kr/kN = 0.8) disordered systems (with 
polydispersity S = 0.01) with F"^^ = Ibrhg {rh is the mean 
particle mass). Thin gray lines correspond to the responses 
of 15 individual realizations in each case, smoothed with a 
Gaussian of width w = 3-R; the thick black line corresponds 
to an average over these 15 configurations. 



lations with polydisperse systems (with radii distributed 
uniformly in the interval [R — S ■ R,R]; see Sec. IIIB 
with S = 0.01,0.1,0.25). In general, we find that poly- 



disperse systems exhibit qualitatively similar behavior to 
the ordered systems described in Sec. |VJ The efi'ects of 
polydispersity are described in detail in this section. 

Since disorder induces fiuctuations (which are not 
present in ordered systems), the forces were coarse 
grained using a Gaussian coarse-graining function in the 
horizontal (x) direction: — L^e-i^/'^) with w = 3R or 

w — 6R, which amounts to calculating the vertical stress 
azz at the floor [9T|. The stress was then averaged over 
several realization of the disorder (typically five) . 

The effect of friction on the response in disordered sys- 
tems is qualitatively similar to that found in the ordered 
lattices described in Sec.jVj Fig. 19 presents a comparison 
of the response (averaged over 15 different realizations, 
coarse grained with w = 3R) obtained in frictionless and 
frictional (^ = 0.2) disordered packings with S = 0.01, 
with an applied force F'^^^ = Ibfhg (rh is the mean par- 
ticle mass). As in the ordered case, the response is very 
different: there are two very distinct peaks in the fric- 
tionless case, but one peak in the frictional case. The 
fluctuations are quite large even for this small degree of 
polydispersity (and larger for the frictionless case), but 
for the magnitude of the applied force used in this case 
the difference between the frictionless and frictional case 
is quite evident even for individual realizations. 



A. Linearity and the Crossover Force 

A notable difference between monodisperse and poly- 
disperse systems is that in the latter, the range of linear 
response is smaller (see Fig. 20 compared to Fig. l6|. 

In order to examine the effects of polydispersity and 
friction in more detail, we performed similar simulations 
with different degrees of polydispersity and coefficients 
of friction /i, for different applied forces. The results are 
depicted in Fig. |2T] The larger the coefficient of friction, 
the larger the applied force at which the crossover from a 
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FIG. 20: The response at a: = (coarse grained with w = 3-R) 
of a single polydisperse systems with 5 — 0.01, normalized by 
the applied force, -Fcxt, vs. the applied force (given in units 
of the mean particle weight, rfig), for different coefficients of 
friction /i, with frictionless walls (/x"''" — 0) and kr/kN = 0.8 
(compare to Fig. |6|. 



single peaked to a double peaked response occurs, much 
like in ordered systems (Sec. VA As mentioned above. 



the crossover force decreases with increasing polydisper- 
sity, or disorder. For polydisperse systems we observed 
a rather gradual crossover even for /i = 1, unlike the 
sharp transition to a highly asymmetric response in lat- 
tice conflgurations [Sec. VE . A possible explanation may 
be that for a single disordered realization, the response 
is not symmetric even for a small external force. In ad- 
dition, it may be impossible to apply very large forces to 
polydisperse systems without causing major rearrange- 
ments. 

It is quite remarkable that an average over a very small 
ensemble of only five realizations (with coarse graining) 
is sufficient to demonstrate the crossover from a single 
peak to a double peak, and even results (in some cases) 
in nearly symmetrical graphs. It appears that, at least 
when the force is not too close to its crossover value, the 
fluctuations within the ensemble are rather small for this 
choice of coarse graining length (see [HI] for a detailed 
study of another set of simulations in the linear force 
regime) , so that the typical response of a realization re- 
sembles that of the ensemble average. 

The results for different degrees of polydispersity (in- 
cluding the case of a lattice) can be summarized in a 
schematic phase diagram. Fig. 22 :F. Note that in 



our simulations a speciflc preparation method has been 
used: relaxation under gravity from a near-lattice con- 
figuration, so that (at least for small 8), the configura- 
tion retains partial order. In the range of polydispersity 
studied here, more disordered systems exhibit a smaller 
crossover force, i.e., they are more susceptible to induced 
anisotropy. It is possible that this result pertains to small 
polydispersity and a larger polydispersity may actually 
stabilize the system. In the systems studied here, the dis- 
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tribution of contact angles may not be quite isotropic, ex- 
hibiting preferred direction close to those of the triangu- 
lar lattice, and therefore more susceptible to failure along 
certain direction. Therefore, in more disordered systems, 
the degree of induced anisotropy may be smaller, yielding 
a non-monotonic dependence of the crossover force on the 
degree of polydispersity. In any case, it is likely that the 
crossover does not depend on the polydispersity alone but 
also on the preparation method (as observed in experi- 
ments on granular piles [15, and slabs [SO])- A more com- 
plete phase diagram would presumably further depend on 
additional parameters which characterize the geometry of 
the packing (in a statistical way). A common approach 
is to use fabric tensors [92], often simplified by consid- 
ering the distribution of contact angles. However, this 
characterization is often employed in conjunction with a 
mean-field approach, which fails for disordered granular 
materials [93j |94] due to the large non-aSine component 
of the microscopic displacements (see also [55]V 



from the point of application as a power law (e.g., the 
solutions for an infinite half plane or half space [Boussi- 
nesq's problem], which decay like l/z in 2D and 1/z^ in 
3D pS', '53^). Therefore, as the distance from the point 
of application increases, the relative displacements of the 
particles become smaller, and they are less likely to reach 
the nonlinear regime (the finite region of open contacts 
discussed in Sec. |V]is an example of such a nonlinear ef- 
fect) , so that far from the point of application of the force, 
the system becomes "more elastic" and the response to 
a distributed load is the same (or nearly the same) as 
the response to the resultant load. This suggestion may 
apply to recent experiments [96, 97 on a 2D system sub- 
ject to a small cyclic displacement, which shows a l/r 
decay of the displacement field at large distances from 
the perturbation. 



VII. INTERPRETATION OF EXPERIMENTAL 
RESULTS 



B. Superposition 

Another issue related to the question of linearity, dis- 
cussed in Sec|VIA| is that of superposition: is the sum of 



responses to several forces applied individually the same 
as the response to all the forces applied at the same time? 

Fig. [23] demonstrates a typical superposition test. 
Here, two downward pointing vertical forces, each of 
magnitude i^oxt = 0.2fhg, are applied at two points at 
the top of the slab, otherwise characterized by S — 1% 
and fi — 0.05, for which the linear range is rather small 
(^cxt ^ OArfig). We compared the response of the system 
when each of the forces was applied separately and simul- 
taneously at different horizontal separations. As shown 
in Fig. |23| superposition holds even on the microscopic 
scale, i.e., for individual forces without coarse graining. 

Interestingly, we find that superposition holds quite 
well even beyond the linear regime, as shown in Fig. [24] 
(for Fcxt = TTig) . A possible reason is that sufficiently far 
from where the forces are applied, the rearrangements 
(changes in the contact network or sliding) are less sensi- 
tive to the precise positions of the applied forces. In more 
general terms, it is possible that the response to a dis- 
tribution of forces, sufficiently far from the region where 
they are applied, is not very sensitive to the details of 
the distribution. 

This explanation is similar to St. Venant's princi- 
ple (see, e.g., |95|). which states (for a linear elastic sys- 
tem) that the difference in stresses and strains in the 
interior of an elastic body due to two separate but stati- 
cally equivalent systems of surface tractions (same over- 
all force and torque), are negligible sufficiently far from 
the area where the loads are applied. In the case con- 
sidered here, the system is certainly not strictly elastic 
for Fcxt ^ QAifig (due to the effects of contact network 
changes and sliding). However, the response to a local- 
ized force (in an elastic medium) decays with the distance 



The results described in the previous section help to 
interpret the experimental results reviewed in Sec. jIIBj 
To this end, it is important to note that the different 
experiments described in Sec. jIIB correspond to differ- 
ent regions of the phase diagram [Fig. 22 described in 
Sec.l^ 

As mentioned, experiments on 2D systems using pho- 
toelastic particles [TS] |22] typically use a rather large 
applied force (i^ext — 150m(7 in these experiments), re- 
quired in order to obtain a significant photoelastic ef- 
fect. In these experiments the interparticle forces, rather 
than the macroscopic stress, were measured (the forces 
on the floor have not been measured). These are repro- 
duced quite well by our DEM simulations with friction 
(see Sec.lv|. It is yet unclear whether these experiments 
are compatible with the phase diagram of Sec. |VI| since 
the macroscopic stress has not been measured directly in 
these experiments. 

The experiments reported in [19l [20] used disordered 
systems with a rather small applied force (a few times the 
particle weight), for which a single peak was observed, 
in agreement with the phase diagram. Since these ex- 
periments used sand, rather than spherical particles, the 
effective coefficient of static friction may be much larger 
than 1, so that sliding is less likely to occur. Note that 
those experiments were performed on a 3D system, while 
the simulations reported here are in 2D; however, we ex- 
pect similar qualitative features of the phase diagram in 
3D systems. The difference in the responses of dense 
and loose sand [50] j as well as the deviations from the 
isotropic elastic prediction observed in both cases, may 
be explained by a (small) anisotropy, induced by the 
packing construction process (see also f55]). We believe 
that the anisotropy induced by the small applied force 
in these experiments should be negligible (note that the 
systems studied in [13] j^H] are also deeper, in terms of 
particle diameters, than the ones studied in [I?ll^ ). 
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FIG. 21: The response of disordered systems with different degrees of polydispersity, 5, and coefficients of friction, (/i = /i"^"), 
with kr/kN = 0.8, to different applied forces T^cxt. The response is averaged over 5 configurations (3 for the largest force used 
with 5 = 10% and /i = 0.2, for which the system was unstable) and smoothed with a Gaussian of width w = 6.R. 




FIG. 22: Schematic phase diagram (in the Fcxt-A* plane) for 
the crossover from a single peaked to a double peaked re- 
sponse, with different degrees of polydispersity. 



In contrast to the experiments reported in [13 1^ . 
the experiments reported in [21] on ordered 3D packings 
used a very large force of a few thousand times the par- 
ticle weights. Several distinct peaks were observed for 
shallow systems. These peaks became more diffuse for 
deeper systems. In the absence of horizontal contacts, 
the coordination number in the HCP and FCC lattice is 
reduced to 6 (each particle is in contact with 3 particles 
in the layer above it and that below it). This corresponds 
to the isostatic limit (in the absence of friction). Appar- 
ently, the experimentalists attempted to avoid horizontal 
contacts by choosing the appropriate wall spacing; even 
if some horizontal contacts were present, it is likely that 
such a large force opened the horizontal contacts, at least 
in the top layers. Since the force is large, it may also lead 
to sliding at most of the contacts (at least in the top lay- 
ers), so that for shallow systems, the observed behavior 
is extremely anisotropic, and, as mentioned, near the iso- 
static limit. However, in deeper systems, the opening of 
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FIG. 23: Superposition in a polydisperse system with S — 0.01, /i = 0.05, and _Fcxt = 0.2fhg, frictionless walls (/i"'^" = 0) and 
kr/kM = 0.8. d indicates the distance between the points of application of the two forces (applied symmetrically with respect 
to the center top particle), given in mean particle diameters R. In the legend, Left denotes the response to a force applied 
to the left of the center particle, Right to a force applied to its right, (Left+Right) to the two forces applied together and 
(Left) + (Right) the sum of Left and Right. The results shown are for the individual forces (no coarse graining) in a single 
realization. 
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FIG. 24: Superposition in a polydisperse system with S = 0.01, /i = 0.05, and Fcxt = rng (same as Fig. 23 with a different 
applied force; the same realization of the disorder was used). 



contacts and sliding may be reduced (as we observed for 
/i — 0.2), reducing the anisotropy (note that even the 
deepest layers in the described experiment were not very 
deep: about 20 layers). In amorphous packings, a single 
peak was observed (however, the method of application 
of the force was different: a force impulse was used, which 
may correspond to a weaker force in our quasi-static de- 
scription). Thus, the results of the experiments reported 
in [21, are also consistent with the schematic phase dia- 



gram presented in Sec. |VI[ 

The phase diagram may also explain the striking dif- 
ference observed in [23' between the experimentally mea- 
sured displacement response in packings of (frictional) 
disks subject to a localized displacement (which exhibits 
a single peak), and the results of simulations of fric- 
tionless isostatic packings (which exhibits two peaks). 
Note, however, that the mean coordination number in 
the frictionless polydisperse systems studied here (e.g., 
z ^ 4.25 for S = 0.01) is well above the isostatic limit 
of z = 4. It is, however, much smaller than that of 
the ideal triangular lattice, z = 6, the "missing" con- 
tacts being predominantly horizontal. This suggests that 
the anisotropic ( "hyperbolic-like" ) response of friction- 
less systems may arise, at least for small polydispersity. 



from the anisotropic structure of these systems rather 
than their (near)-isostaticity, as suggested in [23]. 



VIII. CONCLUDING REMARKS 

We performed simulations studying the response of a 
2D granular slab to a localized force, in particular the 
effects of the magnitude of the applied force, the fric- 
tional parameters and the polydispersity. Our results 
indicate that anisotropy, which is induced by changes in 
the contact network and frictional sliding due to applied 
force, gives rise to a hyperbolic-like response (as already 
suggested in [3l |25]). However, we show that this ef- 
fect becomes less pronounced for large systems (which is 
more likely to be the case in most engineering applica- 
tions) and/or systems subject to smaller applied forces 
for a given system size. This may explain why small- 
scale experiments sometimes yield a hyperbolic-like re- 
sponse while elasticity is typically used for describing 
small deformations in engineering practice. In addition, 
friction, which is present in all real granular systems, in- 
creases the range of linearity in the applied force and 
further reduces the effect of stress-induced anisotropy. 
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rendering the response even closer to that predicted by 
isotropic elasticity. The same behavior is observed in 
polydisperse systems, which is the typical case in nature, 
although in polydisperse systems the forces at which in- 
duced anisotropy becomes large are smaller. The valid- 
ity of superposition beyond the linear elastic regime in- 
dicates that sufficiently far from the applied load, the 
system responds (nearly) elastically even though it is not 
elastic everywhere. 

The above considerations should also be relevant for 
the case of granular piles, whose properties depend 
strongly on their construction history [T5]. For instance, 
when the pile is grown from a point source, it is shaped 
by successive avalanches, all of which are initiated in the 
apex area (beneath the source). The nature of the gran- 
ular flow (i.e., the avalanches) in the case of formation by 
flow from an extended source is clearly different from that 
corresponding to a point source. In particular, the de- 
gree of anisotropy of piles constructed from a point source 
is larger than that of piles built from an extended (and 
approximately uniform) source [99 . In the former case, 
there is a dip in the pressure distribution on the floor, 
whereas in the latter there is no dip. In [TS], this effect 
was described in terms of orientations of force chains, 
which may also reflect the macroscopic anisotropy. How- 
ever, we believe that once the (possibly inhomogeneous) 
anisotropy in the pile is taken into account, the stress 
distribution under the pile should be compatible with an 
elastic picture. The question of how this anisotropy is 
created during construction will probably require a more 
detailed description of the pile dynamics, possibly in- 
volving the changes in the contact network (and sliding) 
which lead to macroscopic anisotropy. 

In general, the elastic description of granular solids is 
expected to be valid for sufficiently large scales (and sys- 
tem sizes) and not-too-large applied forces. Outside this 
range, such a description may need to be extended by us- 
ing a nonlinear, incrementally elastic continuum model, 
with stress-history dependent elastic moduli (which re- 



flect the induced anisotropy, and may exhibit hyperbolic- 
like behavior) . The nonlinear effects of broken contacts 
and frictional sliding should be important in understand- 
ing the failure of granular materials. In particular, the 
dominant effect of tangential forces in such failure, e.g., 
in the transition to an asymmetric response described 
in Se c. |VE| as well as the effects of applied torques 



(Sec. VGl and of rolling resistance, support the sug- 



gestion that, at least on intermediate scales, continuum 
models of micropolar, or Cosserat type [88j [Ml [90] may 
be required (see, e.g., IIOOIIIOIJ I. 

We note, as a suggestion for a future experiment, 
that one should be able to observe experimentally the 
crossover we obtain in the simulations with increasing 
applied force, from a single peaked to a double peaked 
response; in such an experiment, the setup should be 
sufficiently sensitive to measure the response to a small 
force (in the linear regime) but also sufficiently robust to 
enable the application of larger forces. In addition, the 
force would have to be applied slowly, to avoid plastic 
flow of the material. Such an experiment would be very 
useful in testing the applicability of the above physical 
considerations to real granular materials. 
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